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We present a new pseudo-spectral code for the simulation of evolution systems that are second 
order in space. We test this code by evolving a non-linear scalar wave equation. These non-linear 
waves can be stably evolved using very simple constant or radiative boundary conditions, which we 
show to be well-posed in the scalar wave case. The main motivation for this work, however, is to 
evolve black holes for the first time with the BSSN system by means of a spectral method. We 
use our new code to simulate the evolution of a single black hole using all applicable methods that 
are usually employed when the BSSN system is used together with finite differencing methods. In 
particular, we use black hole excision and test standard radiative and also constant outer boundary 
conditions. Furthermore, we study different gauge choices such as 1 -|- log and constant densitized 
lapse. We find that these methods in principle do work also with our spectral method. However, our 
simulations fail after about lOOM due to unstable exponentially growing modes. The reason for this 
failure may be that we evolve the black hole on a full grid without imposing any symmetries. Such 
full grid instabilities have also been observed when finite differencing methods are used to evolve 
excised black holes with the BSSN system. 



PACS numbers: 04.25.Dm, 02.70.Hm, 04.70.Bw, 95.30.Sf 



I. INTRODUCTION 

Currently several gravitational wave detectors such as 
LIGO J,] or GEO are already operating, while several 
others are in the planning or construction phase Q . One 
of the most promising sources for these detectors are the 
inspirals and mergers of binary black holes. In order to 
make predictions about the final phase of such inspirals 
and mergers, fully non- linear numerical simulations of the 
Einstein Equations are required. In order to numerically 
evolve the Einstein equations, at least two ingredients 
are necessary. First we need a specific formulation of 
the evolution equations. And second, a particular nu- 
merical method is needed to implement these equations 
on a computer. For both these ingredients many choices 
are available. However, only a few of these choices have 
so far been successful in the evolution of a binary black 
hole system. In terms of commonly used evolution sys- 
tems there is mainly the BSSN system j3| and also the 
generalized harmonic system The main numerical 
methods are finite differencing and spectral methods. 

The original generalized harmonic system is second or- 
der and has only been used together with finite differ- 
encing techniques IB, H- ^^t, recently a first order 
version of this system has also been successfully evolved 
using a spectral method ^ 0. Nevertheless, the ma- 
jority of the astrophysically relevant binary black hole 
simulations to date have used one particular combina- 
tion: The BSSN evolution system together with finite 
differencing techniques. With this choice, several groups 
have recently performed binary black hole simulations of 
one orbit or more ^Oj IT, 12, 23^ JJ-,, JJi, JJJ- However, 
the BSSN system has never before been tested with a 
spectral method. And in fact, no second order in space 
system such as BSSN has ever been tried together with 
a spectral method. Even for single black holes, only sys- 
tems that are first order in time and space (such as the 



KST system |0| or the first order version of the gen- 
eralized harmonic system) have been used successfully 
with a spectral method. In order to derive such first or- 
der systems additional variables have to be introduced 
which have to satisfy additional constraints. It is there- 
fore possible that second order systems perform better 
since they have fewer potentially unstable constraint vi- 
olating solutions which can be excited due to numerical 
errors. 



The aim of this paper is thus to test the standard BSSN 
system with its standard boundary conditions in a sim- 
ple case but with a pseudo-spectral method. Such a test 
will show how well a spectral method will work when em- 
ployed to this widely used and successful system. We will 
implement the BSSN system in its original second order 
in space form, without introducing any extra variables or 
constraints. This will also allow us to draw some conclu- 
sions about how well second order systems work when a 
spectral method is used. In order to perform this task we 
have developed the SGRID code, which can evolve second 
order in space systems using a pseudo-spectral colloca- 
tion method. This code currently supports cubical or 
spherical domains and allows for spectral expansions in 
Fourier or Chebyshev basis. 



Throughout we will use units where G = c = 1. The 
paper is organized as follows. In Sec. ^ we describe our 
methods and show results for some simple tests where 
we evolve non-linear scalar waves. Sec. IIIII describes our 
particular implementation of BSSN and presents tests for 
case of a single black hole. We conclude with a discussion 
of our results in Sec. IIVI 
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II. DESCRIPTION OF OUR METHOD AND 
CODE TESTS 

In this section we describe the methods we use in 
the SGRID code we have developed. Furthermore, we 
present some code tests with scalar wave equations. 



A. The pseudo-spectral collocation method 

In one spatial dimension, spectral methods are based 
on expansions 



N 



iiX) = ^~aiMX) 



(1) 



1=1 



of every evolved field u{X)m terms of suitable basis func- 
tions Ai{X) with coefficients a;. Once the coefficients are 
known it is easy to compute derivatives of u{X) from 



N 



dxu{x) = ^didxMx), 



(2) 



since the derivatives of the basis functions are known 
analytically. 

However, instead of directly storing and manipulating 

the coefficients a; up to some desired order in Z, we 
make use of the fact that (for the basis functions of in- 
terest) we can derive a; from the values of u{X) at certain 
collocation points. If u{X) is known to have the values 
u{Xi) = Ui at the collocation points X^ for i = 1, 2, N 
it is possible to invert the N equations 
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Y^'^iMXi) 



(3) 
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and to exactly solve for the N coefficients a; in terms of 
the Ui. The location of the different collocation points 
depends on the basis functions used. For example, for 
Fourier expansions the collocation points have to be 
equally spaced in the X-interval considered. 

This approach of storing the field's values Ui at the 
collocation points is called a pseudo-spectral collocation 
method. This method has the advantage that non-linear 
terms such as u{Xi)'^ can be computed from simple mul- 
tiplications. Also, the fields at each point can be evolved 
forward in time by a simple method of lines integrator 
such as Runge-Kutta or iterative Crank-Nicholson (ICN). 

The generalization to 3 dimensions, is straight forward 
and can be summarized by 

u{Xi,Y^,Zk) = dlmnMXi)Bm{Yj)Cn{Zk). (4) 

I.e. we are using basis functions which are products of 
functions that depend only on one coordinate. Note that 
the three coordinates X, Y and Z need not be Cartesian 
coordinates. 



For this, paper we we have chosen X, Y and Z to 
be the spherical coordinates r, 9 and 0. I.e. standard 
Cartesian coordinates are given by 



X = r sm{6) cos{<p) 
y = r sin(^) sin((^) 
z = rcos{0). 



(5) 
(6) 
(7) 



As basis functions we will use the Chebyshev polynomials 

2iVi Rout -^m 



= COS I /arccos 
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(8) 



in the radial direction, where and R.out stand for the 
inner and outer radii of our numerical domain. For the 
angular directions we use the Fourier basis 
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The collocation points are chosen to be 
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(14) 
(15) 
(16) 

Any function u{r, 9, (j)) can then be expressed in this basis 
as 

Nr-l Ne-l N^-l 

u{ri,9j,(j>k) ^ X] X] X] dlmnAl{ri)Bm{9j)Cn{4'k)- 
1=0 m=0 n=0 

(17) 

However, in our code we never really uses this expansion 

to compute all the coefficients di„in- Rather, we only ever 
expand in one direction and instead use 



'■{ri,9j,(f)k) = ^ ai{9j,4>h)Ai{ri), (18) 
;=o 

Ne-l 

'■{ri,9j,(f)k) = ^ bi{ri,^k)Bm{Oj), (19) 

m=0 

i{ri,9j,(f>k) = Yl ^i{ri,9j)Cn{ct>k), (20) 



n=0 



to compute the coefficients di{9j,(j)k), bi{ri,(j)k) or 
Ci{ri,9j) along a line in the r-, 9- or 0-direction. As 
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we will see, this suffices to compute partial derivatives in 
any coordinate direction. 

Note that we always choose Ng to be even, in order 
to assure that collocation points in the 6'-direction (see 
Eq. (|12|l 'l do not occur at = and 9 = tt. In this 
way the coordinate singularities at the poles are avoided. 
Notice that, since both 9 and ip are between and 2tt, we 
have a double covering of the entire domain. This double 
covering is due to the fact that we use Fourier expansions 
in both angles. 

Spherical coordinates also have a coordinate singular- 
ity at r = 0. In this paper we avoid this singularity by 
simply choosing Rm > so that in Eq. is always 
larger than zero. This is the appropriate choice for sim- 
ulating a single black hole using excision. From Eq. ^11^ 
and we see that the collocation points for i = and 
i — Nr~l are located on the inner and outer boundaries 
of the numerical domain. Hence it is straightforward to 
impose boundary conditions there. 



B. Results for non-linear scalar waves 

In order to test our approach we have performed tests 
with a non- linear wave equation for a scalar field ip. 
When written as a first order in time and second or- 
der in space system, the evolution equations in Cartesian 
coordinate are 



grid spacing is of order 



5tV = n 



1 + ^^ 



(21) 
(22) 



where A parametrizes the non-linearity. This is also how 
we have implemented this equation in our computer code, 
with the caveat that the Cartesian derivatives are com- 
puted using the expansions H18|l . (|19|l and H2U|) . By this 
we mean that before each timestep derivatives like dxtp 
are computed using 

where drip, dgip and d^ip can easily be obtained from the 
expansions in Eqs. (|18|l . (|19|l and H20() since the deriva- 
tives of Ai (r) , B„i (9) and C„ (0) are known analytically. 
In order to evolve forward in time we can use any method 
of lines integrator. In this paper we will only present re- 
sults obtained with a second order in time ICN scheme. 
We have, however, also implemented and tested several 
different Runge-Kutta schemes. In terms of stability all 
these integrators work equally well for the cases discussed 
below, so that the ICN results suffice for our purposes. 

The numerical domain for the scalar wave test is a 
spherical shell with an inner radius of — 2 and an 
outer radius of Rout = 52. Since the grid points are 
equally spaced in 9 and (p, the smallest grid spacings 
occur near the poles in the ^-direction. This smallest 



h ~ Rr 



2^ 

'NgN^' 



(24) 



In order to not violate the Currant condition the time 
step At has to be chosen to be of order h as well. This 
means that we have to take very small time steps, so that 
the simulations take longer. But it also means that the 
error in the ICN time integrator will be only of order 

As boundary condition at r = i?j„ we always use 



u{t, Ri, 



= u{t = 0, Rin, 9, < 



(25) 



where u stands for either ip or II. At r — Rout we have 
tested both the constant boundary condition 



u{t, Rout,d, (p) = u{t = 0, Rout, d, 

and also the radiative boundary condition 

u{t, Rout, 6 



(26) 



dtu{t, Rout, Q, < 



-dxu{t, Rout,0, ' 



y 

+ -dyu(t,Rout,0,(j)) 

+ -d,u{t,Rout,e,(P)\.{2r) 

r 1 

For our scalar field example we use v = \, and apply this 
boundary condition to both ip and 11 directly, without 
decomposing ip and 11 into incoming and outgoing modes 
at the boundaries. 

As initial data for the scalar field we use an inward 
moving Gaussian profile given by 



n 
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{x - xq)x [y - yo)y {z - zq)z 



where 

A = 1, = 27, t/o = -zo 0, (Tx 



(28) 



(29) 



10 
(30) 



Figure shows the scalar field energy inside our 
numerical domain, computed from a volume integral over 
the energy density 



II2 
2 



[^2_log(l + V')] 



(31) 



The plot shows for different outer boundary condi- 
tions and also for different A. The dotted line shows the 
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Scalar field energy E 




t 



FIG. 1: The plot shows the scalar field energy inside our 
numerical domain for three different cases. If the fields are 
kept constant at the outer boundary no energy can leave the 
domain and thus the energy remains constant (dotted line). 
If we use the radiative outer boundary condition 12711 . waves 
can leave the domain and thus the energy decreases over time. 
The solid line shows this decrease for the linear (A — 0) case, 
while the broken line depicts the non- linear case (A = 100). 
In all three cases we have used Nr = Ng — 30 and A^'^ — 29. 

result if we use A = 0, together with the outer boundary 
condition (|26|) . In that case the resulting linear waves are 
reflected by the boundaries so that the energy remains 
constant. The solid line depicts the A = case with the 
radiative outer boundary condition H27|) . Since the initial 
wave profile is moving inward, it takes about one cross- 
ing time until the wave starts leaving the grid. After 
that (around t ~ 100) the energy starts dropping. The 
same qualitative behavior occurs in the non-linear case 
for A = 100. Note that even though we used identical 
initial data in all cases, the energy in the non-linear 
case is initially about 35% larger. Thus for A = 100 the 
non-linear terms are not negligible and are by no means 
a small perturbation. This means that the simple radia- 
tive outer boundary condition H27|) is not only capable 
of propagating non-linear waves off the grid, it also does 
not seem to lead to any instabilities in the scalar wave 
tests performed by us. This result is interesting since an 
analogous outer boundary condition is also used in al- 
most all stable black hole evolution codes which use the 
BSSN system together with a finite differencing method. 

In Fig. [21we demonstrate that our code exhibits the ex- 
pected geometric convergence if we increase the number 
of grid points. We show the error in att = 200, for the 
the non-linear wave equation with radiative outer bound- 
ary condition. Since no analytical solution was available 
to us, the error in the energy plotted here is defined by 

Sn,4o{E^) = |i?^,Ar — E^,4o\, (32) 

where E^^pf is the A'^th order spectral approximation to 
E^ obtained by using Nr = Ng = + 1 = N. In Fig. El 



Error in E at t=200 for X =100 
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N 

FIG. 2: This plot shows the error in at t — 200, for the 
non-linear wave equation with the radiative outer boundary 
condition. One can see that the error decreases exponentially 
with the number of grid points N. 



Error in E for ^ =100 with radiative BC 




t 



FIG. 3: This plot shows the error in versus time, for the 
the non-linear wave equation with the radiative outer bound- 
ary condition. 



we show the time evolution of this error. We can clearly 
see how the error decreases with an increasing TV. Figure 
21 shows the scalar field itself at time t = 55 along the 
0-direction, for r = 27 and 9 = n/2. This time was 
chosen such that possible noise form the inner and outer 
boundaries has had time to reach r = 27. Also note that 
during that time the wave which was initially localized 
around (r = 27, 9 = tt/2, (j> — 0) has had time to reach 
the point (r = 27, 9 = 7r/2, (p ~ tt). We see that for 
A^ = 12 points the wave still looks quite jagged, but it 
converges to a smooth result for higher A^. 
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Scalar field \|/(t=55, r=27, 8=71/2, (])) for >.=100 with radiative BC 




1 2 3 4 5 6 



FIG. 4: This plot shows the scalar field if] along an equatorial 
circle {6 — it/ 2) with radius r = 27. Shown are snapshots of 
i/i at t = 55 for different numbers of grid points Nr — Nq — 
N^ + l^N. 

C. Analysis of the well-posedness of the scalar 
wave boundary conditions 

As already mentioned the boundary conditions used 
above are standard, in the sense that they have been 
used in almost all black hole simulations with the BSSN 
system so far. Like in the scalar field example above, they 
are appHed to all evolved fields of the BSSN system, with- 
out decomposing them into incoming or outgoing modes 
at the boundary. Because of this it is not immediately 
obvious whether these boundary conditions are mathe- 
matically well-posed. However, at least in the scalar field 
case well-posedness can be demonstrated easily. 

Scalar wave equations of the form H21|l and H22|l are 
well studied. For the sake of completeness we will here 
quote some results of Gundlach and Martin-Garcia ^3 
who have studied the well-posedness of systems that are 
first order in time and second order in space. First note 
that the non-linearity in Eq. H22|) does not affect well- 
posedness. Hence our Eqs. (PT|) and are completely 
analogous to Eqs. (1) and (2) in According to [l8| . 
the second order characteristic variables for direction rii 
are given by 

U± = n±dnip (33) 
Ua = Oa^, (34) 

where 9„ is the derivative in direction rii and Da are 
the derivatives normal to ni. Here U+ and C/_ are the 
incoming and outgoing modes, while the Ua are zero 
speed modes. By using energy estimates, Gundlach and 
Martin-Garcia ,18] show that the problem is well-posed 
if boundary conditions of the form 

U+ = kU-+ f (35) 

are imposed on the incoming mode only, where |k| < 1 
and / is a given function. 



We will now show that the boundary conditions H26|) 
and H27|l lead to a condition of the form Ij^SI) . 

Let us start with condition (|26l) and assume that 
our initial data satisfy Ii{t — 0, i?out, 6*, (/>) = 
^j{t = 0,Rout,e,(j)) = 0. Then Eq. ^ leads to 
U{t, Rout, 0,(1)) = and t/jit, Rout, 0,(1)) = 0. Now 
Il{t, Rout, 0,(1)) = is completely equivalent to Eq. H35|l 
for K = —1 and / = 0. So if we had only this one 
boundary condition, the problem would certainly be well 
posed. In this case we should compute ip at the boundary 
by using the evolution equation (|21|l . However, if we do 
so, we immediately obtain that ip(t, Rout, 0, (j)) — 0, since 
U{t, Rout, 0,(1)) = 0. Thus imposing ip{t, Rout, 0,(1)) = 
is redundant but completely consistent with the evolu- 
tion equations. Hence our boundary condition (|26|l is 
equivalent to the well posed condition H35|) with k — — 1 
and f — 0. Note that in our actual numerical simu- 
lations we have really set U{t — 0,Rout,0,(j)) — ip{t — 
0,Rout,0,(j)) = 0. I.e. at the boundary we do not really 
have a Gaussian profile. However, since our Gaussian 
is numerically almost indistinguishable from zero at the 
boundary anyway, any discontinuities introduced by this 
procedure are well below our numerical accuracy. 

Next consider the radiative condition (|27|l . For ip it 
implies dtip — —[dni' + 4' I Rout]- If we use only this 
condition for if) together with the evolution equation H21|l. 
we can replace dt'4' by H. In this case wc can compute V' 
from the evolution equation (|21|l and our only boundary 
condition is 

W= - [dn^ + lP/Rout]- (36) 

This boundary condition, however, is again of the form 
(|35|l with K = and / — —ip/Rout- Thus in this case the 
problem would be well posed. 

Now, the second condition coming from Eq. (|27ll . i.e. 

dtii^-[dnii + n/Rout] (37) 

is already implied by the above procedure of using only 
the boundary condition (|36() together with the evolution 
equation H21|l . (To arrive at Eq. H37|l we can simply take 
the time derivative of Eq. H3t)|) and use Eq. (|21|l to replace 
time derivatives of ip.) Hence the radiative conditions 
(|27ll are simply an alternate but equivalent way of setting 
the fields at the boundary. Since this alternate way can 
be derived from the well-posed boundary condition l|3()|l . 
it is also well-posed. 

In summary, the constant as well as the radiative 
boundary conditions of Eqs. iP^ll and l|77|l are both equiv- 
alent to imposing conditions only on the the incoming 
scalar wave modes. Thus the scalar wave system is well- 
posed with either one of the two boundary conditions. 

III. EVOLVING A SINGLE BLACK HOLE WITH 
THE BSSN SYSTEM 

After having demonstrated that our methods lead to 
stable and convergent results for scalar waves we will now 
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turn to evolving a single black hole example using the 
BSSN system ^. 

In the case of BSSN, the 3-metric gij is written as 



1^J 



(38) 



where the conformal metric "fij has unit determinant. In 
addition, the extra variable 



1 z^jJj 



(39) 



is introduced where 7'-' is the inverse of the conformal 
metric. Furthermore, the extrinsic curvature is split into 
its trace free part Aij and its trace K, and given by 



(40) 



In terms of these variables the evolution equation are 



dt<l> 



= -2a A 



£ 



i-aK + A/3' 



-2a ( -^'WjK -QA'WjC, 



(41) 
(42) 



-2A^W,a - e(r - 7^'=f;.,)/3j, + 7-'"'=/3V,z; 



2~ 
3 



=;rp% + p^r^ (43) 



dtA 



dfK 



-i(t> 



-DiDia + a{R 



+ 4) 



TF 



+a{KA,j 

-D'D,a + a ( A'^ A. 



2A,kAf) + £fiA,^ 
1 



(44) 
£pK. (45) 



Here the superscript TF in Eq. denotes the trace 
free part and Di is the derivative operator compatible 
with the 3-metric gij. Notice that F*^ and i?y are the 
Christoffel symbol and Ricci tensor associated with the 
conformal metric 7^ j , while 

Rt = -2b,b,^-2%j&bi(i) 

+Aib,cf>)ibj(t>) ~ A^,,i&cP){Di^), (46) 

where Di is the derivative operator compatible with the 
conformal metric "fij . For all the results reported here we 
have set the constant ^ in Eq. H43() to ^ = 4/3. 

We should also note that in order to ensure that Ay- 
remains traceless during our numerical evolution, we sub- 
tract any trace due to numerical errors after each evolu- 
tion step from Aij. 

As initial data we use the metric and extrinsic curva- 
ture of a Schwarzschild black hole of mass M in Kerr- 
Schild coordinates. Thus initially we have 



9^J 



- S ,. 



(47) 



= a[kHj + IjH^i + Hli^j + HI 



+2H^{lilklj,k + ljhh,k) + "^HliljlkH^k], (48) 



where 

H = M/r 

r = x'/r. (49) 

The indices of V can be lowered with the flat metric. 
As initial lapse and shift we use 



a = l/Vl + 2H, 
13' ^2rH/{l + 2H). 



(50) 
(51) 



If this lapse and shift are kept constant during the evolu- 
tion the right hand sides of Eqs. (|41|) - (|45|l are all zero and 
thus all the BSSN variables are constant. All the results 
presented in this paper are indeed obtained with a con- 
stant shift equal to what is given in Eq. (|51|) . The lapse, 
however, is treated somewhat differentl y a nd we have 
tried two approaches. Firstly, following I20I I21I |2^ 
we have experimented with a so called 1 -I- log-lapse given 
by 



dta = A/3* - aK. 



(52) 



The second approach is to introduce a densitized lapse 
given by 



(53) 



where g is the determinant of the 3-metric. As in [2^ this 
densitized lapse q is then evolved instead of a, which is 
simply computed using Eq. H53() . In our particular case 
we impose 



dtq = 0, 



(54) 



and thus we keep the densitized lapse constant during 
the evolution. Analytically both Eq. and Eq. 
lead to a constant lapse for the Kerr-Schild initial data 
used here. However, as we will see below, numerically 
both cases differ. 

The BSSN system of evolution equations (Eqs. H41|) - 
(|45|l ) is first order in time and second order in space, 
just as the scalar wave example considered before. We 
will thus use the same numerical methods. As our nu- 
merical domain we use a spherical shell, which extends 
from Rin — 1.85Af to Rout — 16A/, with collocation 
points given by Eqs. (|ll|) - (|13|l . On this shell we will use 
the Chebyshev basis ^ in the radial direction and the 
Fourier basis and H10|) in the angular directions. As 
before, we will evolve the Cartesian components of all 
variables on this spherical grid without rewriting the evo- 
lution equations in spherical coordinates. 

The horizon of the black hole is at r = 2M, so that 
the inner boundary at r = Rin — 1.85M is inside the 
horizon, which justifies excising the region r < Rin- We 
treat this inner boundary like a pure outflow boundary 
and do not impose any boundary condition there. Apart 
from using a spectral method this is the one point where 
we really have to do something that differs from what 
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is usually done when the BSSN system is evolved using 
finite differencing methods. Almost all finite differencing 
BSSN simulations performed so far either avoid excision 
altogether or, if excision is used, apply the so called "sim- 
ple excision" boundary condition j2fl], where the time 
derivatives of variables at the inner boundary are ob- 
tained by copying the time derivatives from neighboring 
grid points. Even though one might worry about spoil- 
ing geometric convergence, this "simple excision" algo- 
rithm can in principle also be applied within our spectral 
method using collocation points. The results, however, 
are disastrous and the simulation crashes after only a few 
M of evolution time. In the continuum limit the "sim- 
ple excision" technique corresponds to a Von Neumann 
boundary condition, which can also be applied within our 
spectral approach. This inner boundary condition, how- 
ever, also leads to exponential blow up after only a few M 
of evolution time. Hence the simplest prescription which 
can be used without decomposing our variables into in- 
coming and outgoing modes, is to impose no boundary 
condition at r = Rin- Fortunately this option is easily 
implemented within our spectral approach, and amounts 
to treating the points at r = Rin in the same way as inte- 
rior points by simply using the BSSN evolution equations 
at these points as well. 

At the outer boundary at r = Rout = 16M we have ex- 
perimented with both constant and radiative boundary 
conditions. In the constant case we simply apply Eq. H2ti|) 
to all evolved variables. Analytically this is the correct 
boundary condition since with our choices for lapse and 
shift nothing should evolve, so that all fields indeed re- 
main constant at the boundary. We have also tested 
radiative boundary conditions. Following |2ll | we have 
applied a modified radiative boundary condition to our 
evolved variables. The modification consists of applying 
condition (|27|l to the difference between the evolved vari- 
ables and the analytic solution given by the initial data. 
As in j2l| we have applied this modified radiative bound- 
ary condition to all variables except F*, which was kept 
constant at the outer boundary. 

In the following we will present results for three differ- 
ent combinations of the boundary conditions and lapse 
choices described above. In Fig. [51 we show the ADM 
mass (see e.g. |^) versus time. In principle the ADM 
mass should stay constant for all time. However, due to 
exponentially growing instabilities all our runs eventu- 
ally crash and the ADM mass becomes more and more 
inaccurate over time. The two dashed lines are obtained 
when we use the radiative outer boundary conditions de- 
scribed above. The short dashed line is for the 1 -f log 
lapse of Eq. H52|) . This simulation lasts the longest. The 
long dashed line corresponds to a simulation where the 
densitized lapse of Eq. H53|l is kept constant. The solid 
line is for a constant densitized lapse together with con- 
stant outer boundary conditions. This simulation crashes 
at about the same time as the corresponding one with ra- 
diative boundary conditions (long dashed line). 

Let us discuss next how these results depend on the 
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FIG. 5: ADM mass versus time for different clioices of lapse 
and outer boundary conditions. The short dashed line is for 
the radiative outer boundary condition together with a 1-l-log 
lapse. This simulation crashes a.t t = 137M. The solid line 
depicts the result when all variables are kept constant at the 
outer boundary and a constant densitized lapse is used. In 
this case the run dies at t = 99M. The long dashed line shows 
the case for the radiative boundary condition together with a 
constant densitized lapse. Here a crash occurs at t = IQQM. 
All results shown in this figure were obtained using Nr = 30, 
Ne = 14 and 7V^ = 13. 



Hamiltonian constraint (BC: const, Lapse: densitized, const) 




t/M 



FIG. 6: The plot shows the L'^-norm over the entire grid 
of the normalized Hamiltonian constraint versus time for the 
case of constant outer boundary conditions and a constant 
densitized lapse. In this plot TVe — 1 = A^<^ = 13. For increasing 
Nr we observe geometric convergence in the Hamiltonian until 
t ~ 60 Af. 



resolution used. To this end it is instructive to study 
the normalized Hamiltonian constraint at different reso- 
lutions. Figure El shows the Hamiltonian for an increas- 
ing number of radial grid points Nr in the case of con- 
stant outer boundary conditions and a constant densi- 
tized lapse. After an initial increase in the Hamiltonian 
for the first few M, the curves become flat until insta- 
bilities take over at around t ^ 60 Af. These instabilities 
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Hamiltonian constraint (BC: radiative, Lapse: densitized, const) 
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FIG. 7: The L^-norm of the normahzed Hamiltonian con- 
straint versus time for the case of radiative outer boundary 
conditions and a constant densitized lapse. For increasing A'^,- 
the Hamiltonian converges until about t ~ 60M. The angular 
resolution is iVe — 1 = A^^ = 13. 



eventually lead to a crash. We obtain geometric conver- 
gence until about t ~ 60 Af. Note that the simulations 
last longer when the resolution is increased. 

In Fig. [7|we show the corresponding results for the case 
with radiative outer boundary conditions and a constant 
densitized lapse. Again the result converges with reso- 
lution until about t ^ 60M and simulations run longer 
for higher resolutions. However, this time a flat plateau 
forms only for a short time after the initial increase of 
Hamiltonian. For example for the highest resolution the 
curve is almost flat between 2M and 20M. After that 
the Hamiltonian increases exponentially. This increase is 
the result of constraint violations entering through the 
outer boundary, since the radiative boundary conditions 
used here are not constraint preserving. 

Figure IHl shows our results for the case of 1 + log lapse 
together with radiative outer boundary conditions. The 
qualitative behavior in this case is the same as for the 
constant densitized lapse with radiative outer boundary 
conditions presented before in Fig. |7| Note, however, 
that the time axis has been extended since these runs 
last longer. 

All plots were obtained for an angular resolution of 
Ng = 14 and = 13. All we have done so far was to 
vary the radial resolution. Note that the Schwarzschild 
black hole evolved here is spherically symmetric. Thus, 
when the Cartesian components of our evolved variables 
are expanded in a Fourier series in 9 or (f>, only the three 
lowest frequencies contribute, so that only the first three 
complex coefficients can be non-zero. Since all our vari- 
ables are real, the lowest coefficient must also be real. 
This implies that the first three complex coefficients can 
be described by five real numbers, which corresponds 
to five collocation points. Hence if we assume spheri- 
cal symmetry throughout the evolution, a resolution of 
Ng = = 5 would be entirely sufficient. Recall how- 
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FIG. 8: The L^-norm of the normalized Hamiltonian con- 
straint versus time for the case of radiative outer boundary 
conditions and a 1-l-log lapse. For increasing Nr the Hamilto- 
nian converges until about t ~ 60 Af. The angular resolution 
is iVe - 1 = iV^ = 13. 

ever, that in order to avoid collocation points at the poles 
we require Ng to be even, and thus the minimum we 
need is Ng — 1 ^ = 5. If we run our simulations with 
this minimum angular resolution non-spherical modes are 
suppressed and our simulations run up to 50% longer 
than the ones shown above. Since our ultimate goal is 
to be able to also evolve non-spherically symmetric situ- 
ations we have chosen not to do this. When we increase 
Ng — 1 = N,f, from 5 to a number between 6 and 10 we 
find that run-time decreases. Yet for values of Ng > 14 
and Ncf, > 13 all results are largely independent of an- 
gular resolution. For this reason we have presented only 
results for iVg - 1 = iV^ = 13. 

Recall that in the case of finite differencing BSSN sim- 
ulations with excision, instabilities were observed when 
a single black hole is evolved on a full grid These 
instabilities disappear when symmetries such as octant 
symmetry are imposed. In the present work no sym- 
metries are imposed. This might be one reason for the 
observed instabilities. In fact, here we do not only evolve 
a full grid, but the full grid is also covered twice due 
to our use of Fourier expansions in both angles. This 
double covering corresponds to less symmetries than in 
usual full grid simulations. Thus one might worry that 
it could even lead to additional instabilities, above and 
beyond the full grid instability . This is however not 
the case. 

The double covering can be removed by using the fact 
that all evolved variables u have to satisfy 

u(r, e, 4>) = u{r, 2n~e,(t) + TT) = u{r, 2tt - , (/) - n) (55) 

a TT < 9 < 2tt. Thus if we use an even number of collo- 
cation points in both angles we can simply overwrite all 
variables at collocation points for which tt < 9 < 2TT with 
their values at the corresponding points with < 9 < tt, 
after each time step. Alternatively, one can also take the 
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FIG. 9: The L^-norm of the normahzed Hamiltonian con- 
straint versus time for radiative outer boundary conditions 
and a 1 + log lapse. The thick lines show how the Hamil- 
tonian evolves when we remove the double covering (using 
Eq. 15511 '). For comparison the results with double covering 
(from Fig. |HJ are shown as thin lines. When the double cov- 
ering is removed the runs crash soon after convergence is lost 
(around t ~ 60A/) , while the runs with double covering carry 
on for a while with large errors and without converging. Note 
that in the convergent regime there is no noticeable difference 
between double covering and single covering, and both work 
equally well. The angular resolution is Ne = A^^ — 14. 



average of u at corresponding points and and replace u 
at both points with this average. We have tried both and 
methods to remove the double covering and both give the 
same result. 

The thick lines in Fig. |51 show the time evolution of 
the Hamiltonian constraint for simulations where the the 
double covering has been removed by using the averag- 
ing procedure described above, after each time step. For 
comparison the thin lines also depict the corresponding 
results with double covering. As one can see both sin- 
gle and double covering produce indistinguishable results 
until convergence is lost at around 60M. After that the 
simulations with single covering blow up quickly, while 
the simulations with double covering continue for some 
time (see Fig. (HJ. However, after about 60M the simula- 
tions with double covering have large errors and do not 
converge any longer. Hence the effective runtime during 
which we get useful results is equal for both single and 
double covering. Thus we conclude that both single and 
double covering work equally well and that the double 
covering is not the reason for the observed instabilities in 
the BSSN evolutions described here. 

Another way to remove the double covering would be 
to expand our variables in spherical harmonics instead of 
double Fourier expansions. This is, however, beyond the 
scope of this paper. 

When working with spectral methods it is quite com- 
mon to use filters to suppress instabilities. For example 
in the case of Fourier expansions one often used filter 



consists of zeroing out the top most third of all spectral 
coefficients. This so called 2/3 rule can also be used for 
Chebyshev expansions. We would like to point out that 
no filters were used in the work presented above. We 
have, however, tested the simple 2/3 rule in our case and 
found that it has no effect on the observed instabilities. 
The only effect we observe is a reduced accuracy when 
the filters are on. This reduction is explained by the lower 
effective resolution when 1/3 of all coefficients is zeroed 
out. 



IV. DISCUSSION 

We have, for the first time, evolved the BSSN system 
by means of a pseudo-spectral collocation method. This 
method uses a grid of collocation points that are dis- 
tributed over a spherical shell. On this grid we evolve 
the Cartesian components of all variables. We use Fourier 
expansions in both angles and a Chebyshev expansion in 
the radial direction. This double Fourier expansion leads 
to a double covering of the spherical shell in the sense 
that there are always two grid points which correspond 
to the same physical location in the shell. The fact that 
there are two corresponding points where all evolved vari- 
ables must have the same value can also be used to undo 
the effects of the double covering after each evolution 
step. However, as we have seen in Sec. IHII the simula- 
tions with single and double covering yield basically the 
same results in the convergent regime. 

The BSSN system is second order in space, and we have 
implemented it as is, without any first order reductions. 
The main purpose of this work was to test this system 
together with boundary conditions and gauges that have 
so far only been used together with finite differencing 
methods. In particular, we use black hole excision and 
either radiative or constant outer boundary conditions. 
As gauges we have tested both 1 + log lapse and con- 
stant densitized lapse. Depending on these choices we 
can evolve a single black hole in Kerr-Schild coordinates 
for about lOOM. After that the simulation crashes due to 
unstable exponentially growing modes. While this evo- 
lution time is not too impressive by today's standards, 
it demonstrates that most methods used in codes based 
on finite differencing can in principle also by used with 
spectral methods. In fact, our code tests with a non- 
linear wave equation show that both constant as well 
as radiative boundary conditions can be simply applied 
to all evolved fields without decomposing them into in- 
going and outgoing modes. In this scalar wave case all 
our simulations are stable and no blowups were observed. 
Therefore it seems likely that the observed instabilities 
which occur with the BSSN system when we evolve a sin- 
gle black hole, are not caused by the simple radiative or 
constant outer boundary conditions used here. Rather 
the blowup may due to the use of black hole excision on 
a full grid without imposing any symmetries such as oc- 
tant or bitant symmetry. Such a full grid instability has 
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also been observed in finite differencing implementations 
of BSSN |23| when excision is used. Indeed, we have not 
checked whether the version of the BSSN used here has 
some incoming modes (e.g. gauge modes) at the excision 
boundary. If such modes do exist, they require bound- 
ary conditions at the excision surface and our choice of 
not imposing any boundary conditions there may be the 
cause of the observed instabilities. We leave the study of 
such modes to future work. 

We would also like to point out that the runtime of 
about lOOM in the simulations described above, was 
achieved using the standard BSSN system without any 
parameter fine tuning. For example, it is well known that 
one can modify this system through the addition of con- 
straints on the right hand side of the evolution equations. 
These additional terms can be multiplied by parameters 
which can be adjusted to extend the runtime. Recall that 
when the Einstein-Christoffel system was first imple- 
mented ^3 it ''^^^ fo^' oiily about 40M before it crashed. 
By adding constraints with parameters this system was 
extended into what is now known as the KST system ■ 



After parameter fine tuning this KST system was able 
to evolve single black holes for the much longer time of 
about 8000M [2^. It might be possible to extend the 
runtime of the BSSN simulations described in this work 
in an analogous way. Furthermore, the runtime and ac- 
curacy of our simulations should improve, if the simple 
constant or radiative boundary conditions used here are 
replaced by constraint preserving boundary conditions. 
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